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Abstract: The NNR ( No-Net-Rotation) -MORVEL ( Mid-Ocean Ridge VELocity) 56 is a set 
of angular velocities describing the motions of 56 plates relative to a No-Net-Rotation reference 
frame. These plates can be adjusted in terms of non-overlapping polygonal regions, separated by 
plate boundaries on a unit sphere. During the calculation on the kinematic parameters for these 56 
plates іп a NNR reference frame using the International Terrestrial Reference Frame (ITRF) velocity 
field, the geometric parameters of tectonic plates play a significant role in establishing an absolute 
plate motion model based on space geodesy results. The computational method for these geometric 
parameters implemented as а FORTRAN90 program is described in this paper, allowing an 
evaluation of the area and the inertia tensor of a polygonal region on a unit sphere. This program is 
mainly built on a triangulation algorithm and the adaptive Simpson 's double integral method for 
spherical polygons, which produces highly reliable results for all 56 modern plates. 
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] Introduction 


Most of Earth' s major features can be understood from the interactions between tectonic plates, which 
move independently, separating from, colliding with, and sliding against one another. Until the middle 1960s 
an unifying theory was developed to explain Earth’ s dynamics!" . Several decades after the inception of the 
theory on plate tectonics, the plate dynamic models constructed using the geological and geophysical data have 
been dominant, until long time-span geodetic observations were gathered to estimate contemporary plate 


2531. As one of the most representative geological plate motion models, NUVEL-1A is опе 


kinematic parameters 
of the mainstream models regarding the plate dynamics and kinematics. With the setting up of the enhanced 
amount and quality of the geologic or geodetic data during the last few years, the MORVEL refined the 
precision and accuracy of the geometric and kinematic parameters for 56 plates that are partly taken from an 
updated digital model of pate boundaries by Віті’. Relative to the NUVEL-1A, the MORVEL incorporates 
more than twice as many plates and covers more of Earth' s surface, and nearly all the NUVEL-1A angular 
velocities differ significantly from its MORVEL counterparts" . 

To derive an absolute motion model in a NNR reference frame, however, the inertia tensors are always 
considered as indispensable attribute of all these plates. Despite various established methods for calculating 
plate inertia tensors corresponding to the NUVEL-1A model presented in many papers °” , however it is 
necessary to recalculate а new set for the NNR-MORVELS6 model , given the considerable discrepancy 


between the NUVEL-1A and the MORVEL. 
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When a polygon on the unit sphere is employed for the representation of a tectonic plate，a simplified 
analysis of the plate inertia tensor can be performed through a numerical method, which is carried out over all 56 
plates. The method for calculating all 9 components of the inertia tensor is illustrated in this paper and this 
method requires the precise knowledge of the plate boundaries. The boundary file contains a 2-column sequence 
of the latitude-longitude plate boundary coordinates that fully enclose the plate in the counterclockwise direction. 

In the first section, we introduce some concepts regarding the no-net-rotation conditions and indicate the 
calculation of the Euler vector in an absolute motion model. The following section describes the detailed 
mathematical models to estimate the area and the inertia tensor of the spherical polygons. The last section of 
this paper is dedicated to manifest the results for 56 modern plates and the appendix gives the original 


FORTRAN90 program for obtaining the aforementioned results. 


2 Net Lithosphere Rotation 


А no-net-rotation model for the lithosphere assumes that the integral of v Xr over the Earth’ s surface 
l [u] < 
equals лего, i. e. 


| ox rds =o, (1) 
Earth 


where, r is the radial vector of the surface element on a unit sphere, and v corresponds to the horizontal 


velocity at that position. The angular velocity of net rotation о, was computed as the total angular momentum 


net 


of all plates divided by the moment of inertia of the entire lithosphere, using the equation ^ 9] 


E v хта . (2) 


BT” Earth 


net 


Then it is convenient to convert the Eq.(2) into the following form °! ; 


3 n 
= 5 0o., 3 
о = У, Фо, (3) 


where c; is the Euler vector describing the motion of plate i relative to an inertial reference frame, such as 
ІТКЕ2008, and Q, is the inertia tensor of plate i, where i goes from 1 to n. The angular velocities of the plates 
relative to the NNR reference frame were then found by vector subtraction, namely, 


"i = Ci Е Ow x (4) 


i 


w 
For those tectonic plates where angular velocities are not available in the geodetic model such as the ITRF2000- 
РММ, due to a lack of sufficient data, Altamimi 4! tested four cases to perfect the incomplete geodetic model. 
The fourth case described a method for estimating the missing angular velocity. Here we employ it in this paper 
by using a simple equation written as : 


w +w 


ШЕШ p (5) 
ITRF , А А ds MORVEL . 
where œw; ` is the undetermined rotation vector of plate i in the ITRF model, and c; is the MORVEL 


rotation vector for plate i relative to plate j, which is adjacent to the missing plate i. 


3 Area and inertia tensor of plates 


3.1 Evaluation of the plate area 

The spherical polygons are defined by great circle ares connecting points on the sphere, the positions of 
which are given by latitudes and longitudes. In fact, an algorithm for determining the area of a spherical 
polygon of arbitrary shape has been presented by Bevis and Cambareri  !! , where the kernel idea is to compute 
the interior angle at each vertex of the spherical polygon. In this paper, however, we employ a somewhat 


similar method to that of Miller! trying to determine the area of spherical polygon by summing the signed 
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areas of component triangles. 


For a spherical polygon of n sides, the spherical excess E is generalized as 
E= У a, - ((n- 2) x 180°), (6) 
i=l 


where a,(i=1:…n) are the interior angles of the polygon. Considering a spherical polygon ABCD as shown in 
Fig. 1(a), the north-pole combined with any two adjacent vertices of the polygon can constitute a spherical 
triangle, such as NAB. The two sides of the triangle are known from the latitudes of their vertices, i. e. ап-т/2- 
latitude(A) and bn -m/2-latitude( B). Taking the previously obtained two sides and the included angle specified 
by the difference between longitudes of A and B, the opposite side ab can be calculated via the formula: 


ab = 2sin ! /hav(bn — an) + sinbnsinanhav ZN ， (7) 
where the haversine function is defined as havx = (1-совх)/2 with x in radians. Having obtained all three 


sides of the spherical polygon, we can use the formula to obtain its excess: 


s-an s = bn s— ab 


_ 8 
E = 4tan ' Jian + tan + tan + tan А 
2 2 2 2 


(8) 
an + bn + ab 
2 


To find the area of a spherical polygon, first, one may use the successive vertices in pair to form a 


where s = 


spherical triangle. Each spherical triangle employs the north pole as a common vertex to make the calculations 
convenient. When calculating the areas of the individual triangles, we adopt a convention that the sign of the 
triangle area (which has a same value as the spherical excess for a polygon on a unit sphere) is identical to 
the sign of the difference between the longitudes of a pair of adjacent vertices. If the longitude of the first vertex 
is less than the second one, then the sign of this triangle area is defined as a positive value for a set of points 
arranged in the anticlockwise way and vice versa. Therefore the area of the spherical polygon isregarded as the 
absolute value of the sum of the signed spherical excesses for each of the spherical triangles. Taking the facility 
of the calculation for the upcoming inertia tensor, a provision is crucial that the vertex point traversing the 


polygon prefersto be enumerated in the counterclockwise direction. 


Fig. 1 Spherical polygon and triangles illustrating calculation method for area and inertia tensor discussed in text. 


(a) N triangles constructed from counterclockwise spherical polygon of n sides; 


(b) Spherical triangle encompassing Sorth Pole and the sides of polygon traversing the 180th meridian 


One should note the following two special cases when estimating the geometric parameters of the spherical 


polygons. The first case is as follows; if the sides of the polygon cross over the International Date Line from 
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point E, to W, , as illustrated in the Fig. 1(b) , a 27 has to be added to the longitude of W,. The same applies 
to a similar situation when a vertex leave W, for Е,, a 27 has to be subtracted from E, in order to promise the 
continuity of the calculation. The second case occurs when a spherical polygon has an area larger than that of 
the hemisphere. According to a implicit assumption that the area of any polygon is less Шап 277, for those 
extensive polygons, such as the triangle ABC surrounding the north pole М, the exact area of ABC is the 
complement of ABC including the north pole, i.e. Exactarea = Area‘. = Ат - Area\,.. The original 
FORTRAN90 program in the appendix has taken into account these two singular cases. 

3.2 Estimation of plate inertia tensor 


Тһе components of the symmetric inertia tensor Q сап be calculated for a region P using the following formula: 


0„ = [ pe = x4, (9) 


where x, (4. 7 1,2,3) are the Cartesian coordinates, 0,, (jj, v 2 1,2,3) are the elements of the identity 
matrix, and the integration is carried out over the surface of a plate P. These inertia tensors are based on the 
hypothesis that the surface density of the plate is unit one, and entirely describe the plate geometry. For 


instance, the plate area A is easily calculated by taking the trace of Q: 
100) = У, | Q - xd = 24, (10) 
m 


This implies invariance of the trace under coordinate rotations and the sum of the diagonal components is 
always double the area of the polygon. Generally speaking, non-diagonal components indicate the asymmetry of 
the polygon with respect to the Cartesian axes, and all of the diagonal components have a positive value, which 
is useful, together with the Eq.( 10) , as the verification test for the calculation results. 

In this paper we propose a somewhat different method from Schettino! 7! for constructing the spherical 
triangle. In fact, we will see that the integral at the right-hand side of ( Eq.(9)) is easily calculated for 


spherical triangles. The components of the total tensor are therefore given by: 


О = > | (8, — x,x,) dA 
iz1'm 
| | (11) 
=O 44 = 之 Гах, 
where А is the total polygon area, which has been illustrated in the last section. Let a point on Ше sphere be 


given in spherical coordinates (0,À) , where 0 is the latitude and A is the longitude, so that its Cartesian 


coordinates are given by 
X, = CosOcosA X, = cosOsinA X4 = sin0, (12) 


then the area element dA at this position is equal to cos0dÀd0. The components of the inertia tensor for a 


triangle NAB are therefore written, in spherical coordinates, as: 
А2 7/2 
Q. =б„А-] 4А] Л.Х, 204, (13) 
d 0(A) 


where A, , A, are the longitudes of vertices A, В and the function f is given by 
fu(0, А) fu(0, À) fu(0, А) 
FCO, Л) = |360, А) fa(0, А) fs(0, A) |. (14) 
SaO, А) SaO, А) /ХӨ, А) 
Next, as a result of the symmetry of the inertia tensor, the 6 independent components of f are expressed in the 


following way: 


fa (0, Л) = cos'8cos A (15) 
f (0, A) = cos'0cosAsinA (16) 
5.00, A) = cos20sin0cosÀ (17) 
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fa(0, А) = cos'0sin?A (18) 
fs(0, А) = cos Osin0sinA (19) 
fa (0, А) = sin 0cos0 . (20) 


The upper limit of the inner integral about the latitude is always set to 7/2, because the North Pole is 

considered as the common vertex of each of the spherical triangles. In contrast, with the simple upper limit, 

the lower limit function 0( À ) can be obtained from a serious of derivations, whose concrete form is formulated as 
CicosA + C,sinA 

0(À) =- arctan( с, ) ， (21) 


where C,, C,, C, represent three constants. Once the integrated triangles are determined by one side of the 


polygon, such as AB, we can write their expression in the following form: 
C, = cosÜ,sinA,sinO, - cosO,sinA ,sin6, 
С, = cos0,cosA,sinÜ, - cos0,cosAÀ ,sinO, , (22) 
C, = совб,совб,віп (А, — A,) 

where (0,, A,), (0,, А„) are the latitude and the longitude of vertices A and В, respectively. 

Unlike the process of the area estimation, the evaluation of the inertia tensor is associated with the 
integral order. Hence, a counterclockwise direction must be adopted in the procedure. All of the special 
situations have been considered in the Fortran program , including the 180th meridian case and the encompassing 
the South Pole case. In addition to the aforementioned two special cases in the area estimation, the principal 
moments in the tensor calculation need to be deducted from the whole inertia tensor of the spherical surface in 
order to acquire the exact moments, i.e. Moment;,,, = 8m/3 = Moment principai › when involving the polygon 


containing the south pole. 


4 Results and analysis 


The NNR-NUVEL56 model contains 56 tectonic plates around the earth, and the software OSXStereonet 


developed by Cardozo and Allmendinger' 18: 


was applied to plot the global plate distribution map, as illustrated 
in Fig. 2. Utilizing the Fortran program, we estimated geometry parameters of all 56 modern plates with the 
accuracy the inertia tensor better than 10 ^. Table 1 lists the area and the inertia tensor of all MORVEL plates, 
which provide essential material for calculating plate kinematic parameters in the NNR reference frame. The 
sum of the area of all plates equals 12. 566340 steradian, which is slightly less than the surface area of the 
whole unit sphere, namely 44 (12. 566371) with the relative error 1] = О. 0002396. Our results indicate that the 
relative errors for the six components of the total tensor are 0. 0001796, 0. 0002996 , 0. 0002696, 0. 001096, 
0. 0001696 , and 0. 00010% respectively, It is shown that the inertia tensor of the entire spherical surface is 
8m/3E, where E takes the identity matrix. The discrepancy probably arises from either the imperfection of 


plates over the entire sphere or the unavoidable rounding errors in floating point arithmetic. 


5 Conclusion 


The method for computing the areas and inertial tensors of tectonic plates has been presented. This method 
is based upon the triangulation algorithm and the adaptive Simpson' s double integral procedure, which can be 
applied to the spherical polygons representing such tectonic plates. Results for the NNR-MORVELS6 tectonic 
plates show that highly reliable data can be produced, as long as starting from the precise definition of the plate 
boundaries. In addition, a FORTRAN90 program has been attached to the end of thispaper, which is expected to 


be valuable to the future studies of the kinematics and dynamics associated to the motions of tectonic plates. 
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Fig.2 Plate boundaries and geometries employed for MORVEL 


(a) View direction at 35^ sorth, 0° east; (b) View direction at 20° north, 180° east 


Table 1 Geometric parameters of 56 modern plates included in the NNR-MORVEL56 


Plate" Area" Qu 0 03 Qu Qu Q5 
ram 0. 130659 0. 108248 0. 089481 0. 063589 0. 028732 0. 036320 -0. 051295 
ап 1.432624 1. 326692 1.174711 0. 363845 -0. 050954 0. 052461 0. 081269 
АР 0. 020501 0. 018168 0. 004178 0. 018656 0. 006097 0. 002056 -0. 005420 
аг 0. 120824 0. 074249 0. 066810 0. 100589 -0. 048782 -0. 029553 -0. 031041 
А5 0. 007930 0. 003780 0. 007017 0. 005063 -0. 001938 -0. 003445 -0. 001610 
АТ 0. 014182 0. 008022 0. 011586 0. 008756 -0. 003970 -0. 005767 -0. 003733 
au 0. 921383 0. 597620 0. 558686 0. 686460 0. 223995 -0. 217199 0. 241071 
ВН 0. 012950 0. 007157 0. 005807 0. 012936 0. 006391 -0. 000194 0. 000201 
BR 0. 004814 0. 000353 0. 004786 0. 004488 0. 000322 -0. 001203 0. 000086 
BS 0. 017146 0. 011322 0. 005961 0. 017009 0. 007977 -0. 000850 0. 001173 
BU 0. 012697 0. 012632 0. 000436 0. 012327 0. 000852 0. 000130 -0. 001936 
са 0. 073043 0. 066003 0. 011971 0. 068111 0. 018006 -0. 005213 0. 017059 
СІ. 0. 037650 0. 015349 0. 022487 0. 037464 0. 018192 0. 001581 -0. 001323 
со 0. 072230 0. 071072 0. 003017 0. 070372 -0. 005543 0. 001064 0. 010142 
ер 0. 203647 0. 196537 0. 022175 0. 188580 -0. 021636 0. 007182 0. 045603 
CR 0. 003559 0. 000414 0. 003532 0. 003172 0. 000289 -0. 001100 0. 000101 
EA 0. 004114 0. 003554 0. 001272 0. 003402 -0. 001260 -0. 000631 -0. 001420 
eu 1. 196311 1. 005910 0. 894791 0. 491921 —0. 035559 —0. 213221 -0. 310262 
FT 0. 000789 0. 000054 0. 000787 0. 000736 —0. 000027 -0. 000197 -0. 000007 
СР 0. 000360 0. 000346 0. 000015 0. 000360 -0. 000071 0. 000002 0. 000012 
іп 0. 306360 0. 286350 0. 042318 0. 284052 -0. 057049 -0. 013096 -0. 060490 
jf 0. 006315 0. 005162 0. 004356 0. 003111 —0. 001501 0. 001916 0. 002491 
JZ 0. 002406 0. 002192 0. 000941 0. 001679 -0. 000560 -0. 000394 -0. 001032 
КЕ 0. 012450 0. 003751 0. 012432 0. 008717 -0. 000123 -0. 005589 -0. 000034 
Iw 0. 117084 0. 063137 0. 081094 0. 089937 —0. 043330 0. 036053 0. 029664 
MA 0. 010367 0. 004018 0. 007354 0. 009362 0. 004369 0. 002475 -0. 001708 
MN 0. 000203 0. 000050 0. 000154 0. 000202 0. 000086 —0. 000011 0. 000006 
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Table 1 continued from previous page. 


Plate" Area" Qu Q> Q5 О О Q5 
MO 0. 002841 0. 001271 0. 001581 0. 002830 0. 001405 -0. 000126 0. 000113 
mq 0. 007890 0. 006131 0. 007510 0. 002139 0. 000812 -0. 003172 0. 001465 
MS 0. 010301 0. 007165 0. 003150 0. 010287 0. 004720 -0. 000118 0. 000164 
па 1. 365654 1. 228582 0. 941574 0. 561152 0. 066184 -0. 003632 0. 396247 
NB 0. 009563 0. 002624 0. 006962 0. 009540 0. 004209 -0. 000360 0. 000211 
ND 0. 023942 0. 022362 0. 002000 0. 023523 0. 005755 -0. 000735 0. 002485 
NH 0. 015853 0. 001931 0. 015441 0. 014334 0. 002345 -0. 004547 0. 000758 
NI 0. 003062 0. 000271 0. 003046 0. 002808 -0. 000212 -0. 000839 -0. 000063 
nb 1. 440653 0. 372572 1. 301219 1. 207515 -0. 051346 -0. 005428 0. 044223 
DZ 0. 396683 0. 385354 0. 068410 0. 339603 -0. 012644 -0. 002969 -0. 113428 
OK 0. 074825 0. 053441 0. 066413 0. 029796 0. 012953 0. 030320 -0. 017870 
ON 0. 008000 0. 005617 0. 004075 0. 006307 0. 003044 0. 002004 -0. 002552 
pa 2. 576858 1. 175689 1. 961254 2. 016772 一 0. 429469 0. 077428 -0. 057431 
PM 0. 006744 0. 006575 0. 000333 0. 006580 0. 001000 -0. 000159 0. 001021 
ps 0. 134118 0. 077744 0. 071266 0. 119226 0. 058301 0. 026648 -0. 027639 

ri 0. 002486 0. 002289 0. 000489 0. 002193 —0. 000625 0. 000239 0. 000763 
sa 1. 003382 0. 606780 0. 582701 0. 817282 0. 338318 0. 179187 -0. 168608 
SB 0. 007615 0. 002101 0. 005575 0. 007554 0. 003341 -0. 000571 0. 000346 
SC 0. 041900 0. 036723 0. 034464 0. 012613 0. 005695 0. 011988 -0. 014451 
SL 0. 001780 0. 001675 0. 001490 0. 000396 0. 000174 0. 000381 —0. 000634 
sm 0. 354795 0. 221034 0. 153743 0. 334814 -0. 154899 0. 024755 0. 035861 
sr 0. 027055 0. 018681 0. 026496 0. 008933 0. 001954 0. 012245 -0. 002957 
SS 0. 003170 0. 000715 0. 002505 0. 003119 0. 001275 -0. 000352 0. 000183 
su 0. 219667 0. 188850 0. 036326 0. 214159 0. 069104 0. 006544 -0. 016832 
sw 0. 004543 0. 003525 0. 004269 0. 001292 0. 000527 0. 001817 一 0. 000940 
TI 0. 008704 0. 005784 0. 003120 0. 008503 0. 004009 -0. 000751 0. 001052 
TO 0. 006248 0. 000759 0. 006194 0. 005544 -0. 000536 -0. 001947 -0. 000186 
WL 0. 011163 0. 003492 0. 007835 0. 010998 0. 004966 -0. 001074 0. 000660 
yz 0. 054249 0. 045687 0. 019960 0. 042851 0. 016644 0. 009648 -0. 019528 


“Plate name abbreviations are as follows: am, Amur; an, Antarctic; AP, Altiplano; ar, Arabia; AS, Aegean Sea; au, 
Australia; BH, Birds Head; BR, Balmoral Reef; BS, Banda Sea; BU, Burma; ca, Caribbean; CL, Caroline; cp, Cocos; cp, 
Capricorn; CR, Caroline; EA, Easter; eu, Eurasia; FT, Futuna; GP, Galapagos; in, India; jf, Juan de Fuca; JZ, Juan 
Fernandez; KE, Kermadec; lw, Lwandle; MA, Mariana; MN, Manus; MO, Maoke; mq, Macquarie; MS, Molucca Sea; na, 
North America; NB, North Bismarck; ND, North Andes; NH, New Hebrides; NI, Niuafoou; nb, Nubia; nz, Nazca; OK, 
Okhotsk; ON, Okinawa; pa, Pacific; PM, Panama; ps, Philippine Sea; ri, Rivera; sa, South America; SB, South Bismarck ; 
sc, Scotia; SL, Shetland; sm, Somalia; sr, Sur; SS, Solomon Sea; su, Sundaland; sw, Sandwich; TI, Timor; TO, Tonga; WL, 
Woodlark; yz, Yangtze. Plate abbreviations given in lower case are for plates included in the MORVEL. Plate abbreviations given in 


the upper case are for plates from Bird (2003) ; b. Plate areas are in steradians for a unit sphere, and the sum of which totals 477. 
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The FORTRAN90 Program 
MODULE CALCULATE_SPHPOLAERA_INERTEN 
! 


| THIS IS A FORTRAN90 MODULE, WHICH INCLUDES A SET OF FUNCTIONS AND SUBROUTIBES 

! USED FOR CALCULATING THE AREAS AND INERTIA TENSORS OF AMOUNTS SPHERICAL POLYGONS. 
|жж ck ok ck ck ok ck ck ck ck ck ck ck ck ck ck ck ck ck ck Ck Gk ck ck ck ck ok ok ok ok ck жож ok ok ck Ck ж ok ж 

! DISCUSSION: 

! VERTEXES OF THE SPHERICAL POLYGON MUST BE SEQUENCED IN A COUNTERCLOCKWISE DRIECTION. 
! THE FIRST POINT IN BOUNDARY FILES SHOULD ALWAYS BE IDENTICAL TO THE LAST ONE 

ІІМ ORDER TO CONSTRUCT А CLOSED POLYGON ОҒ N SIDES WITH N+1 DATA POINTS. 

! 

! MODIFIED: 26 DECEMBER 2014 

! 

! AUTHOR; CHUNXIAO LI 


IMPLICIT NONE 
REAL * 8, PARAMETER: : РІ = 3.1415926535897932D0 
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REAL * 8, PARAMETER: : HALFPI = PI/2DO 
REAL * 8, PARAMETER: : DEGREE = 180D0/PI! DEGREES PER RADIAN 
CONTAINS 


! SPITC RETURNS THE AREA AND SIX COMPONENTS OF THE INERTIA TENSOR OF A N-SIDED 
! SPHERICAL POLYGON. 
! THE LAT AND LON SHOULED TAKE DEGREE AS THEIR UNIT, RANGING FROM -90 TO 90 AND 
! FROM -180 TO 180, RESPECTIVELY. 
! THE N+1 DATA POINTS PORTRAY A POLYGON WIYH N VERTEXES. 
| THE OUTPUT PARAMETER, SPHERICALPOLYGONAREA, REPRESENTS THE AREA OF THE SPHERICAL POLYGON 
! IN STERADIANS FOR A UNIT SPHERE. 
! THE REMAMENT SIX PARANETERS INDICATE SIX COMPONENTS OF THE INERTIA TENSOR 
SUBROUTINE SPITC(LAT,LON ,N ,SPHERICALPOLYGONAREA ,SUM11 ,SUM22 ,SUM33 ,SUM12,SUM13,SUM23) 
IMPLICIT NONE 
INTEGER: : N 
REAL * 8:: LAT(N) ,LON(N) 
REAL + 8:: SPHERICALPOLYGONAREA 
REAL = 8:: SUM11,SUM22,SUM33 ,SUM12,SUM13 ,SUM23 
INTEGER: : J 
REAL + 8:: LON1,LON2,LATI , LAT2 
REAL * 8:; HAVB,DLON,PDLON 
REAL * 8: ; T, A, B,C,S, SUM,EXCESS 
REAL * 8;; C1,C2,C3 
COMMON/GROUPI/ C1,C2,C3 
REAL * 8; ; 511,522,533,512,513,523 
REAL + 8: ; EPS 
COMMON/GROUP2/EPS!THE TOLERANCE INVOLVED ON COMPUTATION IS SPECIFIED BY EPS. 


SUM = 0D0 

SUMII = 0D0 

SUM22 = 0D0 

SUM33 = ODO 

SUM12 = ODO 

SUM13 = ODO 

SUM23 = ODO 

DOJ = 1,N-1 
LONI = LON(J)/DEGREE 
LATI = LAT(J)/DEGREE 
LON2 = LON(J+1)/DEGREE 


LAT2 = LAT(J+1)/DEGREE 
CALL COEFFICIENT( LATI ,LON1 ,LAT2 ,LON2) 
PDLON = LON2-LON1 
DLON = ABS(PDLON) 
IF( DLON.GT.1E-6) THEN 
IF( DLON.GT.PI) DLON = 2р0 * PI-DLON 
IF( LON2.LT.LONI.AND.PDLON.LT.-PI) LON2 = LON2 + 2р0 * PI 
IF( LON2.GT.LON1.AND.PDLON.GT.PI) LON2 = LON2 -2D0 * PI 
HAVB = HAV(LAT2-LAT1)+COS(LAT1) * COS(LAT2) + HAV( DLON) 
B = 2D0 * ASIN(SQRT(HAVB) ) 
A - HALFPI-LATI 
C = HALFPI-LAT2 
S = 0.5D0 * ( A+B+C) 
T = TAN(S/2D0) * TAN( (S-A)/2DO) * TAN( (S-B)/2D0) * TAN( (S-C)/2D0) 
EXCESS = ABS(4D0 * ATAN( SORT( ABS(T)))) 
IF ( LON2.LT.LON1) EXCESS = -EXCESS 
SUM = SUM + EXCESS 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT,F11,EPS,S11) 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT, F22, EPS ,S22) 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT , F33 , EPS , S33 ) 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT,F12,EPS,S12) 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT,F13,EPS,S13) 
CALL FSIM2(LONI,LON2,FS LOWERLIMIT,FS UPPERLIMIT , F23 , EPS ,S23 ) 
SUM11 = SUMII + SII 
SUM22 = SUM22 + S22 
SUM33 = SUM33 + 533 
SUMI2 = SUM12 + 512 
SUM13 = SUM13 + 513 
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SUM23 = SUM23 + 523 
END IF 
END DO 
IF(SUM.LT.0.D0) THEN 
SPHERICALPOLYGONAREA = SUM + 4DO * PI 
SUM11 = 8.D0/3.D0 * PI-SUMII-ABS( SUM) 


SUM22 = 8.D0/3.D0 * PI-SUM22-ABS( SUM) 
SUM33 = 8.D0/3.D0 * PI-SUM33-ABS( SUM) 
ELSE 


SPHERICALPOLYGONAREA = SUM 
SUM11 = SPHERICALPOLYGONAREA - SUMII 
SUM22 = SPHERICALPOLYGONAREA - SUM22 
SUM33 = SPHERICALPOLYGONAREA - SUM33 
END IF 
SUM12 = -SUM12 
SUM13 = -SUM13 
SUM23 = -SUM23 
RETURN 
END SUBROUTINE 


! GIVEN THE DOMAIN OF INTEGERAL, FSIM2 RETUENS THE INTEGRAL OF F. 
SUBROUTINE FSIM2(A,B,FS_LOWERLIMIT,FS_UPPERLIMIT,F,EPS,S) 
IMPLICIT NONE 
REAL * 8, EXTERNAL :: FS LOWERLIMIT,FS UPPERLIMIT,F 
REAL +8 :: A,B,EPS,S 
! LOCAL VARIABLES 
REAL «8 :; H,S1,S2, TS1 , TS2, X,G, 80,C 
INTEGER :: N,J 
N=1 
H=0.5D0 * (B-A) 
C=(B-A) * 1.0E-06 
CALL SIMPI(A,FS LOWERLIMIT,FS UPPERLIMIT,F,EPS,S1) 
CALL SIMPI(B,FS LOWERLIMIT,FS UPPERLIMIT,F,EPS,S2) 
Т51-Н ж(51-52) 
DO WHILE ( ABS(H).GE.ABS(C) ) 
Х-А-Н 
TS2=0.5D0 * TS1 
DO J=1,N 
X=X+2D0 * H 
CALL SIMP1(X,FS LOWERLIMIT,FS UPPERLIMIT,F , EPS,G) 
TS2=TS2+H * G 
END DO 
S=(4D0 * TS2-TS1)/3.0D0 
N=N+N 
IF (N.GE.16) THEN 
IF ( ABS(S-S0).LE.EPS * (ABS(S)+1.0D0) ) RETURN 
END IF 
50-5 
TSL=TS2 
H=0.5D0 * Н 
ЕМ” ро 
RETURN 
END SUBROUTINE 


SUBROUTINE SIMP1(X,FS LOWERLIMIT, FS UPPERLIMIT, F, EPS,G) 
IMPLICIT NONE 
REAL * 8 :: X,EPS,G 
REAL * 8, EXTERNAL :: FS LOWERLIMIT, FS UPPERLIMIT, F 
! LOCAL VARIABLES 
REAL * 8 :: Y1, Y2,H, C, TS1 , T2, Y,GO 
INTEGER :: N,I 
N-1 
Y1 = FS LOWERLIMIT( X) 
Y2 - FS UPPERLIMIT( X) 
H=0.5D0 + ( Y2-Y1) 
C- ( Y2-Y1) * 1.0E-06 


天 文 研究 与 技术 35 


TS1=H * (F(X,Y1) HF(X,Y2)) 
DO WHILE (ABS(H).GE.ABS(C)) 
Y zYI1-H 
TS2=0.5D0 * TS1 
DO I- 1,N 
Ү=Ү+2р0 * Н 
TS2- TS24H * F(X,Y) 
END DO 
G=(4D0 * TS2-TS1)/3.0D0 
N=N+N 
IF (N.GE.16) THEN 
IF (ABS(G-GO).LE.EPS * (ABS(G)+1.0D0) ) RETURN 
END IF 
G0=G 
TS1=TS2 
H=0.5D0 * H 
END DO 
RETURN 
END SUBROUTINE 
! FS_LOWERLIMIT RETURN THE LOWERLIMIT CONDUCTING AS A FUNCTION OF 
! LONGITUDE; AND THE FS UPPERLIMIT RETURNS A CONSTANT UPPERLIMIT, РІ/2. 
FUNCTION FS LOWERLIMIT( LON) 
IMPLICIT NONE 
REAL * 8 :: F5 LOWERLIMIT 
REAL *8 :: C1,C2,C3 
COMMON /GROUPI/ СІ, С2,С3 
FS LOWERLIMIT = -ATAN( (Cl * COS(LON) +С2 * SIN( LON) )/C3) 
RETURN 
END FUNCTION 


FUNCTION FS UPPERLIMIT( LON) 
IMPLICIT NONE 
REAL * 8 :: FS UPPERLIMIT 
FS UPPERLIMIT- HALFPI 
RETURN 
END FUNCTION 
! COEFFICIENT RETURNS THREE COEFFICIENTS RELATING THE SPHERICAL COORDINATES OF ANY 
! TWO ADJACENT VERTEXES OF THE POLYGON. 
SUBROUTINE COEFFICIENT( LATI , LON1,LAT2, LON2) 
IMPLICIT NONE 
REAL * 8 :: LATI, LONI,LAT2,LON2 
REAL * 8 :: C1,C2,C3 
COMMON /GROUPI/ СІ, С2,С3 
СІ = COS(LATI) * SIN(LONI) * SIN(LAT2) -COS(LAT2) * SIN(LON2) * SIN(LATI) 
C2 = COS(LAT2) * COS( LON2) * SIN( LAT1)-COS( LATI) * COS(LONI) * SIN(LAT2) 
СЗ = COS(LATI) * COS(LAT2) * SIN(LON2-LON1) 
RETURN 
END SUBROUTINE 


! THE FOLLOWING INTEGRANDS CORRESPOND SIX COMPONENTS OF THE INERTIA TENSOR. 
FUNCTION F11(LON, LAT) 

IMPLICIT NONE 

REAL *8 :: F11 

REAL * 8 :: LAT,LON 

F11 = COS(LAT) * #3 * COS(LON) * *2 

RETURN 
END FUNCTION 


FUNCTION F22( LON, LAT) 
IMPLICIT NONE 
REAL * 8 :: F22 
REAL * 8 :: LAT,LON 
F22 = COS(LAT) ж *3 * SIN(LON) * *2 
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RETURN 
END FUNCTION 


FUNCTION F33(LON,LAT) 
IMPLICIT NONE 
REAL * 8 :: LAT,LON 
F33 = SIN( LAT) ж * 2 ж COS(LAT) 
RETURN 
END FUNCTION 


FUNCTION F12(LON, LAT) 
IMPLICIT NONE 
REAL * 8 :: F12 
REAL * 8 :: LAT,LON 
F12 = COS(LAT) + *3* COS(LON) * SIN(LON) 
RETURN 
END FUNCTION 


FUNCTION F13( LON,LAT) 
IMPLICIT NONE 
REAL * 8 ;; F13 
REAL * 8 :: LAT,LON 
F13 = COS(LAT) * *2* SIN(LAT) * COS(LON) 
RETURN 
END FUNCTION 


FUNCTION F23(LON,LAT) 
IMPLICIT NONE 
REAL * 8 :: F23 
STHZREAL *8 :: LAT,LON 
F23 = COS(LAT) ж * 2 * SIN(LAT) * SIN(LON) 
RETURN 
END FUNCTION 


! HAVERSINE FUNCTION 

FUNCTION HAV(X) 
IMPLICIT NONE 
REAL*8;.X 
REAL *8 :: HAV 
HAV = (1D0-COS( X) )/2D0 
RETURN 

END FUNCTION 
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MORVEL 构造 板块 的 转动 张 量 


1,2 
李 春 晓 
(1. 中 国 科 学 院 云南 天 文 台 ， 云 南 昆明 650011; 2. ФЕЯФЕХФ, ЖЖ 100049) 


摘要 : NNR-MORVEL56 板块 运动 模型 描述 了 全 球 56 个 构造 板块 在 无 整体 旋转 参考 架 下 的 角速度 
运动 参数 。 这 些 板块 可 以 近似 描述 为 单位 球 上 的 无 重 有 登 球面 多 边 形 区 域 。 用 ITRF 速度 场 计算 这 56 个 
板块 相对 于 无 整体 旋转 参考 架 下 的 绝对 运动 时 ， 板 块 的 几何 参数 起 着 至 关 重 要 的 作用 。 详 细 给 出 了 计 
算 板 块 几何 参数 的 方法 并 且 编 写 了 FORTRAN90 程序 以 供 参考 ， 使 得 计算 单位 球 上 板块 的 面积 和 转动 
惯性 张 量 得 以 实现 。 文 中 的 计算 方法 和 程序 主要 采用 球面 三 角 算法 和 自 适应 辛普森 双 积 分 算法 ， 并 对 
AER 56 个 板块 的 几何 参数 进行 了 计算 ， 得 到 了 较为 可 靠 的 计算 结果 。 
关键 词 : 构造 板块 ; 球面 多 边 形 ; 转动 惯量 张 量 ; NNR-MORVELS6 
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